Fusariose: QTLs Description - Part1
Fusariose: QTLs Description - Part1
1/ Introduction
This file aims to describe the QTLs found concerning fusariose resistance.
It concerns the evaluation of every blups computed before.
Let’s upload the file containing all the LOD scores. It contains about 1.6M lines, so it is a bit long
#Watch out, to reproduct analysis, you have to update the path.
data=read.table("~/Dropbox/Publi_Fusariose/ANALYSIS_REPRO/DATA/bilan_simple_marker.gz" , header=T , sep="," )
#data=data[ order(data$LG , data$Distance) , ]
# random lines of the dataset to develop quicker
#data=data[ sample( seq(1,nrow(data)) , 30000 ) , ]
#data=data[ order(data$LG , data$Distance) , ]Let’s upload the genotyping matrix:
# Genotype matrix
geno = read.table(file="~/Dropbox/Publi_Fusariose/ANALYSIS_REPRO/DATA/fichier_genotypage_QTL.csv", sep = ";" , header = F, na.strings = "-")
geno= as.matrix(geno)
colnames(geno)= geno[1,]
geno= as.data.frame(geno[-1 , ])
# Remove useless individuals
#geno <- geno [ geno[,1] %in% names(data),]Let’s upload the phenotyping matrix:
#Watch out, to reproduct analysis, you have to update the path.
pheno=read.table("~/Dropbox/Publi_Fusariose/ANALYSIS_REPRO/DATA/PHENOTYPE/phenotypage_all_fusa_blup.csv" , header=T , sep=";" )
#numeric
pheno[,-1]=apply(pheno[,-1],2,as.numeric)
#check no redondancy in genotype name
table(pheno$geno)[table(pheno$geno)>1]## named integer(0)
# put geno name as rowname
rownames(pheno)=pheno$geno
pheno=pheno[,-1]
# delete columns with only NA
which(apply( pheno , 2 , function(x) all(is.na(x)) )==TRUE)## named integer(0)
pheno=pheno[ , ! apply( pheno , 2 , function(x) all(is.na(x)) ) ]Charge some libraries that will be useful
library(RColorBrewer)
library(rmdformats)
library(xtable)
library(plotly)
library(FactoMineR)
library(knitr)Let’s create a color vector for the study
my_colors = brewer.pal(8, "Set2")
my_colors = colorRampPalette(my_colors)(20)How many variables do we have?
my_carac=levels(data$variable)
nlevels(data$variable)## [1] 115
What is the LOD threshold to declare a QTL significant? I declare it to 3.73
my_QTL_thres=3.73Let’s do some object that contains group of variable:
# reps
REP1=my_carac[grep("rep1", my_carac )]
REP2=my_carac[grep("rep2", my_carac )]
BLUP=my_carac[grep("blup", my_carac )]
# experiments
CAP13=my_carac[grep("CAP13", my_carac )]
GRI11=my_carac[grep("GRI11", my_carac )]
GRI13=my_carac[grep("GRI13", my_carac )]
GRI15=my_carac[grep("GRI15", my_carac )]
LEC14=my_carac[grep("LEC14", my_carac )]
BEQ11=my_carac[grep("BEQ11", my_carac )]
BAR14=my_carac[grep("BAR14", my_carac )]
# Type de phéno NEPI, NEPIL, PEPIL, NOT, DON
PEPIL=my_carac[grep("PEPIL", my_carac )]
PEPI=my_carac[grep("PEPI", my_carac)]
PEPI=setdiff(PEPI, PEPIL)
NEPIL=my_carac[grep("NEPIL", my_carac )]
NOT=my_carac[grep("NOT", my_carac )]
DON=my_carac[grep("DON", my_carac )]
# AUDPC et derniere notation
AUDPC=my_carac[grep("AUDPC", my_carac )]
LAST=c("CAP13_PEPI550_rep1","CAP13_NEPIL550_rep1","CAP13_PEPIL550_rep1","CAP13_PEPI550_rep2","CAP13_NEPIL550_rep2","CAP13_PEPIL550_rep2","CAP13_PEPI550_blup","CAP13_NEPIL550_blup","CAP13_PEPIL550_blup","GRI11_NOT550_rep1","GRI11_NEPIL550_rep1","GRI11_NOT550_rep2","GRI11_NEPIL550_rep2","GRI13_PEPI550_rep1","GRI13_NEPIL550_rep1","GRI13_PEPIL550_rep1","GRI13_NOT550_rep1","GRI13_PEPI550_rep2","GRI13_NEPIL550_rep2","GRI13_PEPIL550_rep2","GRI13_NOT550_rep2","GRI13_PEPI550_blup","GRI13_NEPIL550_blup","GRI13_PEPIL550_blup","GRI13_NOT550_blup","LEC14_PEPI550_rep1","LEC14_NEPIL550_rep1","LEC14_PEPIL550_rep1","LEC14_PEPI550_rep2","LEC14_NEPIL550_rep2","LEC14_PEPIL550_rep2","LEC14_PEPI550_blup","LEC14_NEPIL550_blup","LEC14_PEPIL550_blup","BEQ11_PEPI450","BEQ11_NOT500","BEQ11_NEPIL500","BAR14_PEPI600_rep1","BAR14_NEPIL600_rep1","BAR14_PEPIL600_rep1","BAR14_PEPI600_rep2","BAR14_NEPIL600_rep2","BAR14_PEPIL600_rep2","BAR14_PEPI600_blup","BAR14_NEPIL600_blup","BAR14_PEPIL600_blup")
# agronomical variable
EPI=my_carac[grep("_EPI" , my_carac )]
FLO=my_carac[grep("FLO" , my_carac )]
HEI=my_carac[grep("HEI" , my_carac )]
VER=my_carac[grep("VER" , my_carac )]
# Resistance set for publication = BLUps when avail. DO we have all the needed traits?
a=c(DON, NOT, PEPI, NEPIL, PEPIL)
b=intersect(a,BLUP)
c=intersect(BEQ11, a)
PUBLI=c(b,c,"LEC14_DON")Let’s check what is available for which experiment?
# start table
bilan=data.frame( matrix(0,5,7) )
colnames(bilan)=c("CAP13", "GRI11", "GRI13", "GRI15", "LEC14", "BEQ11","BAR14")
rownames(bilan)=c("DON?", "PEPI?", "NOT?", "NEPIL?", "PEPIL?")
all=list(CAP13, GRI11, GRI13, GRI15, LEC14, BEQ11, BAR14)
num1=0
num2=0
for(i in all){
num1=num1+1
num2=0
for(j in list(DON, PEPI, NOT, NEPIL, PEPIL)){
num2=num2+1
a=intersect(i, j)
bilan[ num2 , num1]=length(a)
}}
kable(xtable(bilan))| CAP13 | GRI11 | GRI13 | GRI15 | LEC14 | BEQ11 | BAR14 | |
|---|---|---|---|---|---|---|---|
| DON? | 1 | 0 | 1 | 0 | 0 | 0 | 0 |
| PEPI? | 4 | 0 | 4 | 3 | 4 | 3 | 5 |
| NOT? | 0 | 6 | 4 | 3 | 0 | 5 | 1 |
| NEPIL? | 4 | 6 | 4 | 3 | 4 | 5 | 6 |
| PEPIL? | 4 | 0 | 4 | 3 | 4 | 0 | 5 |
We know that:
- DON is available for CAP13, GRI13 and BAR14 only. - PEPI is not available for GRI11 - NOT is not available for CAP13. - PEPIL is not available for GRI11 and BEQ11. So please check it is indeed what we have in this QTL file!
2/ Useful functions
A set of function allowing to characterize significant QTLs for a set of variable:
#Same function but with interactive graphics!
plot_a_interactiveQTL=function(data, select_chromo, my_variable, LOD_threshold, my_ylim){
data=data
select_chromo="4B"
my_variable=HEI
p=data %>%
filter(LG==select_chromo & variable %in% my_variable) %>%
group_by(variable, Distance) %>%
select(Distance, variable, LOD, R2, group_physique, Posi_physique) %>%
summarise(
LOD=max(LOD),
R2=max(R2, na.rm=T),
Posi_physique=mean(Posi_physique, na.rm=T),
group_physique=table(group_physique)[1]
) %>%
mutate(my_text=paste("<br>", "Pos phy: ", group_physique , " | ", Posi_physique , sep="")) %>%
ggplot(aes(x=Distance, y=LOD, color=variable)) +
geom_line()
# interactive
ggplotly(p)
}
#A function to plot a QTL and visualize the LOD scores along the chromosome. Works for several variable on a specific chromosome.
plot_a_QTL=function(data, select_chromo, variable, LOD_threshold, my_ylim){
# order data
data=data[ order(data$LG , data$Distance) , ]
# plot the first variable
a=data[data$variable==variable[1] & data$LG==select_chromo & !is.na(data$LOD) , ]
plot(a$LOD ~ a$Distance , type="l" , ylim=c(0,my_ylim) , ylab="LOD score" , xlab="position (cM)" , col=my_colors[1] , lwd=1.6 )
# for every other variable, add a LOD score
num=1
for( i in c(2:length(variable))){
num=num+1
a=data[data$variable==variable[num] & data$LG==select_chromo & !is.na(data$LOD) , ]
points(a$LOD ~ a$Distance , type="l" , ylim=c(0,my_ylim) , ylab="LOD score" , xaxt="n" , col=my_colors[num] , lwd=1.6 )
}
# add lod threshold
abline(h=my_QTL_thres , col="grey")
# add legend
legend("topright" , horiz=F , col=my_colors[c(1:length(variable))] , legend=variable , bty="n" , lty=1 , pt.cex=2, lwd=1.6)
legend("topleft" , horiz=F , col="white" , legend=select_chromo , bty="n" , lty=1 , pt.cex=0, lwd=1.6)
}
#A function to summarize all the significant QTL of a list of variables. Gives the max LOD score, the IC, the position etc..
summary_table_for_QTL=function(data, selected_variable, LOD_threshold, size_IC){
# order data
data=data[ order(data$LG , data$Distance) , ]
my_colors=brewer.pal(5,"Paired")
# create empty summary table
bil=data.frame(matrix(0,1000,10))
colnames(bil)=c("pop","carac","chromo","LOD_max","position","marker","IC","R2","a","diff")
# loop to study every variable & chromosome
num_line_bilan=0
num_chromo=0
for(chrom in levels(data$LG)){
num_variable=0
num_col=0
num_chromo=num_chromo+1
for(var in selected_variable){
num_col=num_col+1
current_data=data[ data$variable==var & data$LG==chrom , ]
current_data=current_data[!is.na(current_data$LOD) , ]
signif_data=current_data[ current_data$LOD > LOD_threshold , ]
#Si j'ai des marqueurs significatifs
if(nrow(signif_data)>0){
#Je récupère les infos de ce QTL
LOD_max=max(signif_data$LOD)
mark_max=signif_data$marqueur[signif_data$LOD==LOD_max]
pos_max=signif_data$Distance[signif_data$LOD==LOD_max]
r2_max=signif_data$R2[signif_data$LOD==LOD_max]
a_max=signif_data$a[signif_data$LOD==LOD_max]
my_diff=ifelse( signif_data$moy.A[signif_data$LOD==LOD_max][1] > signif_data$moy.B[signif_data$LOD==LOD_max][1] , "Dic2>Silur" , "Silur>Dic2" )
#Détermination IC --> je bloque la zone max a 30 cM
in_IC=current_data[ current_data$LOD > (LOD_max-size_IC) & current_data$Distance > (pos_max-30) & current_data$Distance < (pos_max+30) , ]
IC_min=min(round(in_IC$Distance,2))
IC_max=max(round(in_IC$Distance,2))
#Je remplie le tableau bilan
num_line_bilan=num_line_bilan+1
bil[num_line_bilan , 1]="DS"
bil[num_line_bilan , 2]=var
bil[num_line_bilan , 3]=chrom
bil[num_line_bilan , 4]=LOD_max
bil[num_line_bilan , 5]=pos_max[1]
bil[num_line_bilan , 6]=as.character(mark_max)[1]
bil[num_line_bilan , 7]=paste(IC_min, IC_max, sep="-")
bil[num_line_bilan , 8]=r2_max[1]
bil[num_line_bilan , 9]=a_max[1]
bil[num_line_bilan , 10]=my_diff
#close 3 loops
}}}
# Clean and print the Table
bil=bil[bil$pop!=0 , ]
bil=bil[ order(bil$pop , bil$chromo , bil$carac) , ]
return(bil)
#close function
}A set of function allowing to study confidence intervals:
#A function to represent confidence interval of QTLs on the genetic map.
show_IC_QTL_on_map=function(data, selected_variable, selected_chromosome, LOD_threshold, size_IC){
# user can pick every variable
if(selected_variable[1]=="all"){selected_variable=levels(data$variable)}
# order data
data=data[ order(data$LG , data$Distance) , ]
my_colors=brewer.pal(12,"Set3")
my_colors = colorRampPalette(my_colors)(30)
var_with_QTL=c()
# loop to study every variable & chromosome
num_chromo=0
# prepare background of chart
nchro=length(selected_chromosome)
xchro=c(seq(5,5*nchro,5) )
par(mar=c(0,5,2,1))
plot(1,1,col="transparent",bty="n", xaxt="n", yaxt="l" , xlab="", ylab="position (cM)" , xlim=c(3,max(xchro)+11) , ylim=c(320,-5) )
num=0
for(chrom in selected_chromosome){
num=num+1
A=data[data$variable==selected_variable[1] & data$LG==chrom, ]
points(rep(xchro[num],nrow(A)) , A$Distance , pch=20 , cex=0.8 , col="grey")
}
# IC calculation and plot
for(chrom in selected_chromosome){
num_variable=0
num_chromo=num_chromo+1
num_col=0
for(var in selected_variable){
num_col=num_col+1
current_data=data[ data$variable==var & data$LG==chrom , ]
current_data=current_data[!is.na(current_data$LOD) , ]
signif_data=current_data[current_data$LOD>LOD_threshold , ]
#Si j'ai des marqueurs significatifs
if(nrow(signif_data)>0){
# I add it to the list of variables with QTL
if(!var%in%var_with_QTL){var_with_QTL=c(var_with_QTL,var)}
# Je récupère les infos de ce QTL
LOD_max=max(signif_data$LOD)
mark_max=signif_data$marqueur[signif_data$LOD==LOD_max]
pos_max=signif_data$Distance[signif_data$LOD==LOD_max]
# Détermination IC --> je bloque la zone max a 30 cM
in_IC=current_data[ current_data$LOD > (LOD_max-size_IC) & current_data$Distance > (pos_max-30) & current_data$Distance < (pos_max+30) , ]
IC_min=min(in_IC$Distance)
IC_max=max(in_IC$Distance)
# J'ajoute le trait corresponant a ma variable
num_variable=num_variable+0.2
lines(c(xchro[num_chromo]+num_variable,xchro[num_chromo]+num_variable) , c(IC_min,IC_max) , col=my_colors[which(var_with_QTL%in%var)], lwd=6)
num_variable=abs(num_variable)
#Close 2 loops and a if
}}}
#Ajout légende et nom de chromosome?
text(xchro , rep(-10,3) , selected_chromosome , col="orange")
#legend("bottomright", 14 , 240 , horiz=F , col=my_colors[1:length(var_with_QTL)] , legend=var_with_QTL , bty="n" , lty=1 , lwd=6 )
#Close function
}
#A function that gives IC features: size, # of markers, # of genes...
find_IC_infos=function(min_IC, max_IC, chrom, to_remove){
# Get the part of the genetic map concerned by this IC:
IC_map=data[ which(data$LG==chrom & data$Distance>=min_IC & data$Distance<=max_IC & data$LG==data$group_physique & data$variable==levels(data$variable)[1]), c(1:5)]
# remove markers if needed.
if(length(to_remove)>0){
IC_map=IC_map[-which(IC_map$marqueur%in%to_remove) , ]
}
# Calculate statistics
# genet
deb_IC_gen=min(IC_map$Distance)
end_IC_gen=max(IC_map$Distance)
size_IC_gen=end_IC_gen - deb_IC_gen
size_IC_gen=round(size_IC_gen , 2)
nb_of_mark=nrow(IC_map)
# physic
deb_IC_phy=min(IC_map$Posi_physique)
end_IC_phy=max(IC_map$Posi_physique)
size_IC_phy=end_IC_phy - deb_IC_phy
nb_of_gene=nrow(POS[which(POS$group==chrom & (POS$position*1000000)>=deb_IC_phy & (POS$position*1000000)<=end_IC_phy) , ])
# return a vec with all info
vec=c(chrom, min_IC, max_IC, size_IC_gen, nb_of_mark, deb_IC_phy, end_IC_phy, size_IC_phy, nb_of_gene)
return(vec)
}
#A function to show all the genes present in a IC. The relatioship between physical and genetic position is also available
show_genes_of_IC=function(min_IC, max_IC, chrom, my_lwd, my_cex, to_remove){
# get data for the IC
IC_map=data[ which(data$LG==chrom & data$Distance>=min_IC & data$Distance<=max_IC & data$group==data$group_physique & data$variable==levels(data$variable)[1] & data$LG==data$group_physique ) , c(1:5)]
# remove markers if needed.
if(length(to_remove)>0){
IC_map=IC_map[-which(IC_map$marqueur%in%to_remove) , ]
}
# physical informations
#deb_IC_phy=min(IC_map$Posi_physique)
#end_IC_phy=max(IC_map$Posi_physique)
my_genes_pos=POS[which(POS$chromo_BW==chrom) , ]
#plot
plot(IC_map$Distance ~ as.numeric(IC_map$Posi_physique/1000000),
ylim=c(min_IC, max_IC),
xlab="physical pos. of markers (Mb)", ylab="genetic position (cM)" ,
main=paste("chromome ",chrom ," -- IC : ",min_IC," - ",max_IC," cM",sep="" ),
col=rgb(0.4,0.3,0.6,0.4) , pch=20 , cex=my_cex
)
abline(v=my_genes_pos[,3]/1000000 , col="grey", lwd=my_lwd)
}What happens on a specific position
# A function to represent a QTL. It makes a boxplot with 2 categories: alleles A and B
boxplot_a_QTL=function(marker, trait, ...){
# prepare data: merge expression and alleles:
nucl= geno[,c(1,which(colnames(geno)==marker)) ]
my_pheno= data.frame(ind=rownames(pheno), pheno[ , trait])
AA=merge(my_pheno, nucl, by.x=1, by.y=1, all=T)
colnames(AA)=c("Genotype", "Phenotype", "Allele")
AA$Allele[which(AA$Genotype=="dic2")]="A"
AA$Allele[which(AA$Genotype=="silur")]="B"
AA=na.omit(AA)
AA$Genotype
# Highlight Dic2 and Silur
AA$type=factor(ifelse(AA$Genotype=="dic2" | AA$Genotype=="silur","Highlighted",AA$Allele))
# Make plot
p=ggplot(AA, aes(x=Allele, y=Phenotype, text=Genotype, fill=type, color=type )) +
geom_boxplot(alpha=0.5) +
geom_jitter() +
theme(legend.position="none") +
ggtitle(trait)
# interactive
ggplotly(p)
}What about epistasic relationships
#A function to calculate interaction between 2 QTLs for a trait
# Function to calculate interaction effect between QTLs:
analyse_inter=function(marker1, marker2, trait){
# prepare data: merge expression and alleles:
nucl= geno[ , c(1, which(colnames(geno)%in%c(marker1,marker2)) ) ]
my_pheno= data.frame(ind=rownames(pheno), pheno[ , trait])
AA=merge(my_pheno, nucl, by.x=1, by.y=1, all=T)
AA=na.omit(AA)
# Prepare
my_mark=paste(AA[,3],AA[,4],sep="-")
means <- round(tapply(AA[,2],my_mark,mean,na.rm=T) ,2)
# Calculate interaction
model=lm(AA[,2] ~ AA[,3] * AA[,4])
sum=summary(model)
tot_r2=round(sum$r.squared,3)
inter_pval=round(sum$coefficients[4,4],10)
# Complete summary files
res=c(trait, marker1, marker2, means, tot_r2, inter_pval)
return(res)
}
#A function to boxplot 2 QTLs together
# A function to represent a QTL. It makes a boxplot with 2 categories: alleles A and B
boxplot_two_QTL=function(marker1, marker2, trait, ...){
# prepare data: merge expression and alleles:
nucl= geno[ , c(1, which(colnames(geno)%in%c(marker1,marker2)) ) ]
my_pheno= data.frame(ind=rownames(pheno), pheno[ , trait])
AA=merge(my_pheno, nucl, by.x=1, by.y=1, all=T)
AA=na.omit(AA)
# Prepare
my_mark=paste(AA[,3],AA[,4],sep="-")
means <- round(tapply(AA[,2],my_mark,mean,na.rm=T) ,2)
#plot
par( mar=c(5,5,2,2))
boxplot(AA[,2] ~ my_mark ,
medlwd=0, cex.axis=0.6, cex.col="grey", las=1, main="" ,
...,
ylab=trait, xlab="Bi-locus genotype" , col=my_colors[c(3,6,8,3)] , xaxt="n", boxwex=0.4)
my_labels=c( expression(paste('Dic2'["1B"],"-",'Dic2'["5A"])) ,
expression(paste('Dic2'["1B"],"-",'Silur'["5A"])) ,
expression(paste('Silur'["1B"],"-",'Dic2'["5A"])) ,
expression(paste('Silur'["1B"],"-",'Silur'["5A"])) )
axis(labels=my_labels , at=c(1,2,3,4) , side=T)
}
# A function to represent a QTL. It makes a boxplot with 2 categories: alleles A and B
boxplot_two_QTL_interactive=function(marker1, marker2, trait, ...){
# example to develop if needed
#marker1="Cluster_16778|Contig1|original@1840"
#marker2="Cluster_9940|Contig1|complementarySeq@104"
#trait="BAR14_PEPIL300_blup"
# prepare data: merge expression and alleles + add dic2 and Silur
nucl= geno[ , c(1, which(colnames(geno)%in%c(marker1,marker2)) ) ]
qq=data.frame(c("dic2","silur"), c("A","B"), c("A","B"))
colnames(qq)=colnames(nucl)
nucl=rbind(nucl,qq)
my_pheno= data.frame(ind=rownames(pheno), pheno=pheno[ , trait])
AA=merge(my_pheno, nucl, by.x=1, by.y=1, all=T)
AA=na.omit(AA)
# Prepare
AA$my_mark=as.factor(paste(AA[,3],AA[,4],sep="-"))
means <- round(tapply(AA[,2],AA$my_mark,mean,na.rm=T) ,2)
AA=AA[order(AA$my_mark),]
AA$text=paste("ind: ",AA$ind,"<br>","value: ",round(AA$pheno,2),sep="")
# First we calculate a new x axis with the jiiter added to it. Jitter level must be relative to the number of points per level.
my_prop=summary(AA$my_mark) / nrow(AA)
new_x=c()
for(i in c(1:nlevels(AA$my_mark))){
myjitter<-jitter( rep(i, summary(AA$my_mark)[i]), amount=my_prop[i]/2)
new_x=c(new_x, myjitter)
}
# Plot
plot_ly(AA) %>%
add_trace( x=as.numeric(AA$my_mark), y = ~pheno, color = ~my_mark, type = "box" ) %>%
add_trace( x=new_x, y = ~AA$pheno,
type = "scatter",
mode="markers",
marker = list(size = 10, color = ifelse(AA$ind=="dic2","blue",ifelse(AA$ind=="silur","red",rgb(0.1, 0.1, 0.1, 0.3))), line = list(color = rgb(0.1, 0.1, 0.1, 0.8), width = 2)),
text=~text,
hoverinfo="text"
)%>%
layout(title = "",
hovermode="closest",
showlegend=F,
xaxis = list(title = "" , tickfont = list( color="orange", size=26 ), tickmode="array", tickvals=c(1,2,3,4) , ticktext=levels(AA$my_mark) , showline=T ),
yaxis = list (title = colnames(data)[2] , gridwidth=2, zeroline=F)
)
}3/ Determine LOD threshold
We used R-QTL to do some permutation and determine which is the LOD threshold over which q QTL can be considered as significant. (This is done on a cluster because of long computation time). A lod threshold was calculated for each trait. The highest treshold observed was 3.73. We thus used this threshold for every trait since it is the most conservative one.
cd ~/work/FUSA/LOD_SEUIL
cd /NAS/g2pop/HOLTZ_YAN_DATA/DIC2_SILUR/QTL/PUBLI
# je créé un script avec:
for i in $(seq 2 318) ; do echo $i ; Rscript /NAS/g2pop/HOLTZ_YAN_DATA/programmes/Calculate_LOD_threshold.R genotypage.csv phenotypage.csv carte $i ; done
# puis j'envoie
qsub -q normal.q -b yes -cwd -N LOD_thres "./script"
more LOD_thres.o* | grep -B1 -A2 "1000 permutations)" | grep -v "1000 perm" | grep -v "lod" | sed 's/\[1\] \"//' | sed 's/\"//' | sed 's/5% //' > LOD_threshold_per_trait
cat LOD_threshold_per_trait | grep -v "-" | grep -v "[A-Z]" > tmp
R
data=read.table("tmp")
max(data[,1])
min(data[,1])
hist(data[,1])
4/ Agronomical features
Let’s check if we found the expected QTLs for height and flowering date (already known in the litterature..).
height
Supposed to be on chromosome 4B. As you can see on the graphic below, there is no doubt concerning this QTL. It also proove that scripts are working properly. Note you can zoom and hover the graphics to get more details!
plot_a_interactiveQTL(data, "4B", HEI , my_QTL_thres, 40)Same result in a table:
summary_QTL_HEI=summary_table_for_QTL(data, HEI, my_QTL_thres, 1.5)
kable(xtable(summary_QTL_HEI))| pop | carac | chromo | LOD_max | position | marker | IC | R2 | a | diff |
|---|---|---|---|---|---|---|---|---|---|
| DS | BEQ11_HEI | 4B | 10.07464 | 57 | Cluster_3278|Contig2|original@188 | 57-57 | 11.33570 | 6.912055 | Dic2>Silur |
| DS | CAP13_HEI_blup | 4B | 21.42926 | 57 | Cluster_3278|Contig2|original@188 | 57-57 | 36.53591 | 0.000000 | Dic2>Silur |
| DS | GRI13_HEI_blup | 4B | 67.63065 | 57 | TRIDC4BG006730.1@987 | 57-57 | 64.86658 | 0.000000 | Dic2>Silur |
| DS | GRI15_HEI_blup | 4B | 11.73699 | 57 | Cluster_3278|Contig2|original@188 | 57-57 | 24.19335 | 0.000000 | Silur>Dic2 |
We can see that this QTL is really consistent from one experiment to another. Dic2 is taller than Silur. Let’s visualize it in a boxplot:
boxplot_a_QTL("TRIDC4BG006730.1@987", "CAP13_HEI_blup" )flowering (heading?) date
Flowering date and precocity are known to be controlled by a QTL on chromosome 7B. We also find this result:
p1=plot_a_interactiveQTL(data, "7B", EPI , my_QTL_thres, 30)## Warning: We recommend that you use the dev version of ggplot2 with `ggplotly()`
## Install it with: `devtools::install_github('hadley/ggplot2')`
p2=plot_a_interactiveQTL(data, "7B", FLO , my_QTL_thres, 30)## Warning: We recommend that you use the dev version of ggplot2 with `ggplotly()`
## Install it with: `devtools::install_github('hadley/ggplot2')`
subplot(p1, p2)Verse
There is only one QTL for the “verse”, on the same position that the QTL that controls height. This is quite logical: the higher is the plant, the more chance it has to fall… Verse is available for CAP13 only. It is funny to see that the QTL is significant for rep1 and rep2, but not for the blup…
plot_a_interactiveQTL(data, "4B", VER , my_QTL_thres, 30)## Warning: We recommend that you use the dev version of ggplot2 with `ggplotly()`
## Install it with: `devtools::install_github('hadley/ggplot2')`
Epaisseur des glumes
TODO
5/ Resistance to fusariose
Toxine concentration (DON)
Let’s summarize in a table every significant QTL for the “DON” variable:
summary_QTL_DON=summary_table_for_QTL(data, DON, my_QTL_thres, 1.5)
kable(xtable(summary_QTL_DON))| pop | carac | chromo | LOD_max | position | marker | IC | R2 | a | diff |
|---|---|---|---|---|---|---|---|---|---|
| DS | CAP13_DON_blup | 3A | 3.930642 | 168.5 | TRIDC3AG066720.1@1104 | 164.4-175.3 | 10.4574 | 95.60727 | Silur>Dic2 |
If we consider only blups, there is only one QTL on chromosome 3A, with a LOD of 4.19 for DON in CAP13.
Let’s show this 3A chromosome with the DON variables:
plot_a_interactiveQTL(data, "3A", intersect(DON,PUBLI), my_QTL_thres, 6)## Warning: We recommend that you use the dev version of ggplot2 with `ggplotly()`
## Install it with: `devtools::install_github('hadley/ggplot2')`
There is absolutely no effect for other trials than CAP13…
Type 1: Note FUSA
Summary of all the significant QTLs for the ** visual notation**:
summary_QTL_NOT=summary_table_for_QTL(data, NOT, my_QTL_thres, 1.5)
kable(xtable(summary_QTL_NOT))| pop | carac | chromo | LOD_max | position | marker | IC | R2 | a | diff |
|---|---|---|---|---|---|---|---|---|---|
| DS | BAR14_NOT600_blup | 1B | 4.976054 | 1.8 | TRIDC1BG001890.1@27 | 0-13.2 | 11.5233980 | 0.1354486 | Dic2>Silur |
| DS | BEQ11_NOT450 | 1B | 6.309231 | 0.5 | TRIDC1BG001340.4@2565 | 0.5-5.1 | 14.5990816 | 0.6443723 | Dic2>Silur |
| DS | BEQ11_NOT500 | 1B | 7.822870 | 0.5 | TRIDC1BG001340.4@2373 | 0.5-1.8 | 17.7999965 | 0.6567100 | Dic2>Silur |
| DS | BEQ11_NOTAUDPC | 1B | 4.893092 | 0.5 | Cluster_1882|Contig7|original@596 | 0-6.2 | 0.9858993 | 106.7617271 | Dic2>Silur |
| DS | GRI11_NOT450_blup | 1B | 3.936324 | 142.9 | Cluster_16753|Contig2|original@619 | 136.8-146.4 | 10.2581917 | 0.1430546 | Silur>Dic2 |
| DS | GRI13_NOT550_blup | 1B | 4.254277 | 160.4 | TRIDC1BG068250.1@1950 | 160.4-173.9 | 9.1647534 | 0.1740732 | Silur>Dic2 |
| DS | GRI15_NOT350_blup | 1B | 6.450241 | 5.1 | TRIDC0UG001570.3@1699 | 1.3-14.5 | 14.4900190 | 0.0673185 | Dic2>Silur |
| DS | GRI15_NOT550_blup | 1B | 7.803936 | 1.3 | TRIDC1BG000060.1@264 | 0.5-6.2 | 17.3751226 | 0.2101685 | Dic2>Silur |
| DS | GRI15_NOTAUDPC_blup | 1B | 8.064012 | 5.1 | TRIDC0UG001570.3@2193 | 1.3-6.2 | 17.8945425 | 43.2173606 | Dic2>Silur |
| DS | GRI13_NOT350_blup | 3A | 4.777256 | 183.7 | TRIDC3AG069520.1@1240 | 175.3-195.5 | 10.3227644 | 0.0691291 | Silur>Dic2 |
| DS | GRI13_NOT450_blup | 3A | 4.139495 | 183.7 | TRIDC3AG069520.1@672 | 173.7-190 | 8.9076580 | 0.0777861 | Silur>Dic2 |
| DS | GRI13_NOTAUDPC_blup | 3A | 5.142387 | 183.7 | TRIDC3AG069520.1@672 | 175.3-190 | 11.1175468 | 53.7728680 | Silur>Dic2 |
| DS | BEQ11_NOT300 | 5A | 6.468229 | 239.3 | Cluster_265|Contig1|likelySeq@300 | 236.9-239.3 | 11.3610513 | 0.3236671 | Silur>Dic2 |
| DS | BEQ11_NOT350 | 5A | 7.519043 | 237.5 | Cluster_149|Contig3|original@533 | 233.9-239.3 | 17.1751744 | 0.5422521 | Silur>Dic2 |
| DS | BEQ11_NOT450 | 5A | 5.981729 | 247.2 | TRIDC5AG069790.7@31 | 237.5-248.3 | 13.8769478 | 0.6393258 | Silur>Dic2 |
| DS | BEQ11_NOTAUDPC | 5A | 7.708366 | 237.5 | TRIDC5AG067550.4@2885 | 236.9-247.2 | 1.5192131 | 146.0311422 | Silur>Dic2 |
| DS | GRI15_NOT550_blup | 7B | 4.645063 | 166.2 | TRIDC7BG065740.3@1878 | 157.1-166.2 | 10.5634901 | 0.1409060 | Dic2>Silur |
I have 2 main QTLs: on the 1B and on the 5A.
Then I have some really small messages on 3A for Gri13, and 7B for Gri15
Let’s show the QTLs on graphics. We show only the last note and audpc for blup when available.
toshow=c("GRI11_NOT550_blup","GRI11_NOTAUDPC_blup","GRI13_NOT550_blup","GRI13_NOTAUDPC_blup","GRI15_NOT550_blup","GRI15_NOTAUDPC_blup","BEQ11_NOT500","BEQ11_NOTAUDPC","BAR14_NOT600_blup")
#Chromosome 1B
plot_a_interactiveQTL(data, "1B", intersect(NOT,PUBLI), my_QTL_thres, 10)## Warning: We recommend that you use the dev version of ggplot2 with `ggplotly()`
## Install it with: `devtools::install_github('hadley/ggplot2')`
#Chromosome 5A
plot_a_interactiveQTL(data, "5A", intersect(NOT,PUBLI), my_QTL_thres, 10)## Warning: We recommend that you use the dev version of ggplot2 with `ggplotly()`
## Install it with: `devtools::install_github('hadley/ggplot2')`
Let’s visualize an example trough a boxplot:
a=boxplot_a_QTL("TRIDC0UG001570.3@2193", "GRI15_NOTAUDPC_blup" )## Warning: We recommend that you use the dev version of ggplot2 with `ggplotly()`
## Install it with: `devtools::install_github('hadley/ggplot2')`
b=boxplot_a_QTL("TRIDC5AG067550.4@2885", "BEQ11_NOTAUDPC")## Warning: We recommend that you use the dev version of ggplot2 with `ggplotly()`
## Install it with: `devtools::install_github('hadley/ggplot2')`
subplot(a,b)Type 1: % Spike with Fusa
A lot of significant QTLs.
summary_QTL_PEPI=summary_table_for_QTL(data, PEPI, my_QTL_thres, 1.5)
kable(xtable(summary_QTL_PEPI))| pop | carac | chromo | LOD_max | position | marker | IC | R2 | a | diff |
|---|---|---|---|---|---|---|---|---|---|
| DS | BAR14_PEPI350_blup | 1B | 7.383115 | 5.1 | TRIDC0UG001570.3@2193 | 0.5-5.1 | 16.798434 | 7.1975506 | Dic2>Silur |
| DS | BAR14_PEPI400_blup | 1B | 6.446925 | 5.1 | TRIDC0UG001570.3@1806 | 0.5-14.5 | 14.813339 | 6.9793838 | Dic2>Silur |
| DS | BAR14_PEPI550_blup | 1B | 8.231304 | 11.9 | TRIDC1BG003330.11@3771 | 0.5-13.2 | 18.524351 | 10.6648096 | Dic2>Silur |
| DS | BAR14_PEPI600_blup | 1B | 11.514558 | 1.8 | TRIDC1BG001560.1@330 | 1.3-1.8 | 24.614519 | 2.9335351 | Dic2>Silur |
| DS | BAR14_PEPIAUDPC_blup | 1B | 9.760591 | 1.3 | TRIDC1BG001970.2@434 | 1.3-6.2 | 21.473043 | 3897.5655867 | Dic2>Silur |
| DS | GRI13_PEPI450_blup | 2B | 4.436633 | 14.6 | TRIDC2BG002330.3@1130 | 12.2-18 | 9.571098 | 0.2884493 | Dic2>Silur |
| DS | GRI13_PEPIAUDPC_blup | 2B | 3.788355 | 14.6 | TRIDC2BG002330.3@1130 | 12.2-18 | 8.114487 | 506.7832162 | Dic2>Silur |
| DS | GRI13_PEPI350_blup | 3A | 5.065082 | 183.7 | TRIDC3AG069520.1@672 | 180.1-184.9 | 10.950229 | 1.9672612 | Silur>Dic2 |
| DS | GRI13_PEPI450_blup | 3A | 3.899620 | 184.9 | Cluster_1251|Contig2|complementarySeq@569 | 175.3-190 | 8.366828 | 0.2872906 | Silur>Dic2 |
| DS | GRI13_PEPIAUDPC_blup | 3A | 4.980949 | 183.7 | TRIDC3AG069520.1@1240 | 179-190 | 10.767507 | 610.7678883 | Silur>Dic2 |
| DS | BEQ11_PEPI350 | 4A | 4.074615 | 107.6 | Cluster_948|Contig2|original@2756 | 105.9-107.6 | 8.783936 | 6.4687500 | Dic2>Silur |
| DS | BEQ11_PEPIAUDPC | 4A | 4.030530 | 107.6 | TRIDC4AG045960.1@749 | 105.9-107.6 | 8.613270 | 3208.9948529 | Dic2>Silur |
| DS | BEQ11_PEPI350 | 4B | 4.112312 | 48.8 | TRIDC4BG005390.2@1159 | 47.6-69.2 | 8.868261 | 7.0819805 | Silur>Dic2 |
| DS | BAR14_PEPI350_blup | 5A | 5.890239 | 229.1 | TRIDC5AG065190.1@501 | 229.1-237.5 | 13.593976 | 5.9203939 | Silur>Dic2 |
| DS | BAR14_PEPI400_blup | 5A | 5.044484 | 229.1 | TRIDC5AG065260.2@750 | 229.1-238 | 11.680486 | 6.7378080 | Silur>Dic2 |
| DS | BAR14_PEPI600_blup | 5A | 4.693709 | 237.5 | TRIDC5AG067550.4@2885 | 233.9-239.3 | 10.867117 | 1.8876568 | Silur>Dic2 |
| DS | BAR14_PEPIAUDPC_blup | 5A | 6.422906 | 229.1 | TRIDC5AG065260.2@750 | 229.1-237.5 | 14.762259 | 2854.9994215 | Silur>Dic2 |
| DS | BEQ11_PEPI350 | 5A | 5.310014 | 237.5 | Cluster_343|Contig2|likelySeq@420 | 236.9-257.6 | 11.480355 | 8.7632477 | Silur>Dic2 |
| DS | BEQ11_PEPIAUDPC | 5A | 5.086288 | 237.5 | Cluster_149|Contig5|original@1002 | 233.9-257.6 | 10.911406 | 4406.6703131 | Silur>Dic2 |
| DS | GRI13_PEPI550_blup | 5A | 3.916479 | 252.4 | TRIDC5AG071590.5@815 | 246-277.5 | 8.404146 | 1.7975920 | Silur>Dic2 |
| DS | GRI15_PEPI550_blup | 5A | 4.707260 | 242.1 | Cluster_11094|Contig1|complementarySeq@313 | 232.1-245.6 | 10.644227 | 3.7144168 | Silur>Dic2 |
| DS | GRI15_PEPIAUDPC_blup | 5A | 4.639353 | 242.1 | Cluster_13017|Contig1|original@433 | 231.3-245.6 | 10.488543 | 781.0038900 | Silur>Dic2 |
| DS | CAP13_PEPIAUDPC_blup | 6A | 3.929697 | 157.5 | TRIDC6AG056520.1@2331 | 157.5-163.3 | 10.311553 | 826.7490335 | Silur>Dic2 |
| DS | BEQ11_PEPI450 | 7B | 4.395495 | 39.0 | TRIDC7BG009820.4@918 | 16.8-47.2 | 9.192045 | 8.5826325 | Dic2>Silur |
| DS | CAP13_PEPI400_blup | 7B | 8.080811 | 16.8 | TRIDC7BG001540.1@113 | 16.8-16.8 | 17.034043 | 11.3752637 | Silur>Dic2 |
| DS | GRI15_PEPI550_blup | 7B | 4.684315 | 157.9 | Cluster_6415|Contig2|original@249 | 157.1-166.2 | 10.591689 | 3.4147478 | Dic2>Silur |
| DS | GRI15_PEPIAUDPC_blup | 7B | 3.975488 | 166.2 | TRIDC7BG065740.3@489 | 157.1-168.4 | 8.944878 | 622.8496833 | Dic2>Silur |
Two main QTLs on chromosome 1B and 5A. For the 1B QTLs, it is interesting to see that once more, the resistance allele comes from Silur.. How is it possible? The 1B QTLs works for BAR14 only, but for every traits and with strong LOD-scores.The 5A works for 3 experiments.
Then we have some small signals on chromosomes 2A, 2B, 3A, 4A, 6A. Nothing really interesting.
And once more the QTL linked with height and precocity on chromosomes 4B and 7B. For the 4B–> no BLUP.
We can have a look to the LOD-scores of the main traits (AUDPC and LAST):
# Variable to show on the plot?
my_var_to_show=c( Reduce(intersect, list(PEPI,BLUP,AUDPC)), Reduce(intersect, list(PEPI,BLUP,LAST)) )
par(mfrow=c(2,2))
#Chromosome 1B
plot_a_interactiveQTL(data, "1B", my_var_to_show, my_QTL_thres, 12)## Warning: We recommend that you use the dev version of ggplot2 with `ggplotly()`
## Install it with: `devtools::install_github('hadley/ggplot2')`
#Chromosome 5A
plot_a_interactiveQTL(data, "5A", my_var_to_show, my_QTL_thres, 7)## Warning: We recommend that you use the dev version of ggplot2 with `ggplotly()`
## Install it with: `devtools::install_github('hadley/ggplot2')`
An example for QTL 1B and 5A on a boxplot:
Let’s visualize an example trough a boxplot:
a=boxplot_a_QTL("TRIDC1BG001970.2@434", "BAR14_PEPIAUDPC_blup")
b=boxplot_a_QTL("TRIDC5AG065260.2@750", "BAR14_PEPIAUDPC_blup")
subplot(a,b)Type 2: % of Spikelet with Fusa
A lot of significant QTLs.
summary_QTL_PEPIL=summary_table_for_QTL(data, PEPIL, my_QTL_thres, 1.5)
kable(xtable(summary_QTL_PEPIL))| pop | carac | chromo | LOD_max | position | marker | IC | R2 | a | diff |
|---|---|---|---|---|---|---|---|---|---|
| DS | BAR14_PEPIL300_blup | 1B | 5.082209 | 1.3 | TRIDC1BG001970.2@352 | 0.5-6.2 | 11.768093 | 0.9102612 | Dic2>Silur |
| DS | BAR14_PEPIL400_blup | 1B | 5.209982 | 1.8 | TRIDC0UG002680.1@311 | 0-14.5 | 12.061141 | 1.6308273 | Dic2>Silur |
| DS | BAR14_PEPIL550_blup | 1B | 10.983280 | 1.3 | TRIDC1BG000060.1@264 | 0.5-5.1 | 23.688770 | 5.1209416 | Dic2>Silur |
| DS | BAR14_PEPIL600_blup | 1B | 47.428998 | 1.3 | TRIDC1BG001970.2@582 | 1.3-1.3 | 58.793684 | 20.7970099 | Dic2>Silur |
| DS | BAR14_PEPILAUDPC_blup | 1B | 15.252221 | 1.3 | TRIDC1BG001970.2@502 | 1.3-1.8 | 30.560263 | 1499.2994809 | Dic2>Silur |
| DS | CAP13_PEPIL300_blup | 3B | 3.869784 | 128.6 | Cluster_312|Contig3|original@253 | 108.7-131.6 | 8.907305 | 0.9340312 | Silur>Dic2 |
| DS | BAR14_PEPIL300_blup | 5A | 3.845809 | 237.5 | Cluster_343|Contig2|likelySeq@420 | 219.6-247.2 | 8.850201 | 0.8289052 | Silur>Dic2 |
| DS | BAR14_PEPIL550_blup | 5A | 3.915950 | 229.1 | TRIDC5AG065260.2@750 | 220.6-246 | 9.019606 | 2.8962405 | Silur>Dic2 |
| DS | CAP13_PEPIL400_blup | 7B | 5.319279 | 16.8 | TRIDC7BG001540.1@113 | 16.8-18.7 | 12.866356 | 4.8895228 | Silur>Dic2 |
QTL 1B is present for BAR14 only, once more with the Silur allele beeing resistant. The 5A QTL is once more available, with the Dic2 allele beeing resistant. Small QTL on 3B. 4B and 7B QTL are present as usual but nor for BLUPs!
We can have a look to the LOD-scores:
# Variable to show on the plot?
my_var_to_show=c( Reduce(intersect, list(PEPIL,BLUP,AUDPC)), Reduce(intersect, list(PEPIL,BLUP,LAST)) )
#Chromosome 1B
plot_a_interactiveQTL(data, "1B", my_var_to_show, my_QTL_thres, 45)## Warning: We recommend that you use the dev version of ggplot2 with `ggplotly()`
## Install it with: `devtools::install_github('hadley/ggplot2')`
#Chromosome 5A
plot_a_interactiveQTL(data, "5A", my_var_to_show, my_QTL_thres, 7)## Warning: We recommend that you use the dev version of ggplot2 with `ggplotly()`
## Install it with: `devtools::install_github('hadley/ggplot2')`
Type 2: Nb of Spikelet with Fusa
A lot of significant QTLs.
summary_QTL_NEPIL=summary_table_for_QTL(data, NEPIL, my_QTL_thres, 1.5)
kable(xtable(summary_QTL_NEPIL))| pop | carac | chromo | LOD_max | position | marker | IC | R2 | a | diff |
|---|---|---|---|---|---|---|---|---|---|
| DS | BAR14_NEPIL300_blup | 1B | 4.450629 | 1.8 | TRIDC0UG002680.1@315 | 0-11.9 | 10.295939 | 4.8652348 | Dic2>Silur |
| DS | BAR14_NEPIL400_blup | 1B | 4.751407 | 1.3 | TRIDC1BG001970.2@458 | 0.5-14.5 | 11.001903 | 8.6271445 | Dic2>Silur |
| DS | BAR14_NEPIL550_blup | 1B | 10.625259 | 1.3 | TRIDC1BG000060.1@264 | 0.5-5.1 | 23.052405 | 26.6781814 | Dic2>Silur |
| DS | BAR14_NEPIL600_blup | 1B | 47.399229 | 1.3 | TRIDC1BG000060.1@264 | 1.3-1.3 | 58.778490 | 102.8526378 | Dic2>Silur |
| DS | BAR14_NEPILAUDPC_blup | 1B | 13.706392 | 1.8 | TRIDC0UG002680.1@315 | 1.3-1.8 | 28.215512 | 7778.0452855 | Dic2>Silur |
| DS | BEQ11_NEPIL450 | 1B | 9.653514 | 1.3 | TRIDC1BG001970.2@596 | 0.5-1.8 | 5.980977 | 37.0416667 | Dic2>Silur |
| DS | BEQ11_NEPIL500 | 1B | 7.761044 | 1.3 | TRIDC1BG001970.2@457 | 0.5-1.8 | 2.257181 | 34.1111111 | Dic2>Silur |
| DS | BEQ11_NEPILAUDPC | 1B | 7.964935 | 1.8 | TRIDC1BG001970.2@423 | 0.5-5.1 | 4.195497 | 4895.1983447 | Dic2>Silur |
| DS | GRI13_NEPIL550_blup | 1B | 3.791953 | 173.9 | TRIDC1BG070800.4@550 | 160.4-189.4 | 8.122491 | 9.3823938 | Silur>Dic2 |
| DS | GRI11_NEPIL550_blup | 2B | 4.144849 | 108.5 | Cluster_18588|Contig1|original@301 | 98.6-110 | 10.816610 | 3.7106381 | Dic2>Silur |
| DS | BEQ11_NEPIL300 | 5A | 8.214926 | 237.5 | Cluster_149|Contig3|original@533 | 237.5-239.3 | 17.186178 | 6.5306102 | Silur>Dic2 |
| DS | BEQ11_NEPIL350 | 5A | 7.324755 | 237.5 | Cluster_149|Contig5|original@1002 | 233.9-239.3 | 8.609625 | 18.0060217 | Silur>Dic2 |
| DS | BEQ11_NEPIL450 | 5A | 4.687992 | 245.6 | Cluster_6432|Contig1|original@938 | 237.5-250.1 | 3.053695 | 28.2070833 | Silur>Dic2 |
| DS | BEQ11_NEPILAUDPC | 5A | 6.987653 | 247.2 | TRIDC5AG069790.7@31 | 237.5-248.3 | 3.725141 | 4648.8848315 | Silur>Dic2 |
| DS | GRI11_NEPIL450_blup | 5A | 4.005286 | 46.0 | TRIDC5AG005440.2@840 | 39.7-47.6 | 10.443455 | 0.6073008 | Silur>Dic2 |
| DS | GRI11_NEPIL500_blup | 5A | 3.854576 | 47.6 | Cluster_7428|Contig3|original@3241 | 36.4-53.4 | 10.037839 | 1.1677623 | Silur>Dic2 |
| DS | BAR14_NEPIL_blup | 7A | 5.877454 | 185.3 | Cluster_2718|Contig1|original@181 | 185.3-208.3 | 13.563564 | 12.0307361 | Dic2>Silur |
| DS | CAP13_NEPIL400_blup | 7B | 5.054352 | 16.8 | TRIDC7BG001540.1@113 | 16.8-18.7 | 12.235941 | 8.8491561 | Silur>Dic2 |
The 1B QTL is present for BAR14 as usual, with Silur as resistant. Warning, GRI13 as a QTL on the 1B also, but not on the same position.
Small QTL on 1A, 2B, 3A, 3B, 4A, 6B, 7A. QTL height and preco as usual, but not for BLUP.
We can have a look to the LOD-scores:
# Variable to show on the plot?
my_var_to_show=c( Reduce(intersect, list(NEPIL,BLUP,AUDPC)), Reduce(intersect, list(NEPIL,BLUP,LAST)) )
par(mfrow=c(2,2))
#Chromosome 1B
plot_a_interactiveQTL(data, "1B", my_var_to_show, my_QTL_thres, 45)## Warning: We recommend that you use the dev version of ggplot2 with `ggplotly()`
## Install it with: `devtools::install_github('hadley/ggplot2')`
#Chromosome 5A
plot_a_interactiveQTL(data, "5A", my_var_to_show, my_QTL_thres, 7)## Warning: We recommend that you use the dev version of ggplot2 with `ggplotly()`
## Install it with: `devtools::install_github('hadley/ggplot2')`
6/ Summary table
Per experiment
Let’s summarize what QTL has been found for each experiment. For each expe and each phenotype, we write QTL found, and the number of phenotypes that are significant for this QTL between bracklets.
# initialize table
bil=data.frame(matrix(0,5,7))
colnames(bil)=c("CAP13", "GRI11", "GRI13", "GRI15", "LEC14", "BEQ11", "BAR14")
rownames(bil)=c("DON", "NOT","PEPI", "NEPIL", "PEPIL")
# fill it
# For each variable (DON, EPIL,...)
num=0
for(j in list(summary_QTL_DON, summary_QTL_NOT, summary_QTL_PEPI, summary_QTL_NEPIL, summary_QTL_PEPIL)){
num=num+1
# And for each experiment
num2=0
for(i in colnames(bil)){
num2=num2+1
a=j[ grep(i, j$carac ) , ]
b=table(a$chromo)
c=paste(names(b), " (", b, ")", sep="")
d=paste(c, collapse = " - ")
bil[num,num2]=d
}}
# write when data is not available
bil[1,c(2,4,6,7)]="no data"
bil[2,c(1)]="no data"
bil[3,c(2)]="no data"
bil[5,c(2,6)]="no data"| CAP13 | GRI11 | GRI13 | GRI15 | LEC14 | BEQ11 | BAR14 | |
|---|---|---|---|---|---|---|---|
| DON | 3A (1) | no data | () | no data | () | no data | no data |
| NOT | no data | 1B (1) | 1B (1) - 3A (3) | 1B (3) - 7B (1) | () | 1B (3) - 5A (4) | 1B (1) |
| PEPI | 6A (1) - 7B (1) | no data | 2B (2) - 3A (3) - 5A (1) | 5A (2) - 7B (2) | () | 4A (2) - 4B (1) - 5A (2) - 7B (1) | 1B (5) - 5A (4) |
| NEPIL | 7B (1) | 2B (1) - 5A (2) | 1B (1) | () | () | 1B (3) - 5A (4) | 1B (5) - 7A (1) |
| PEPIL | 3B (1) - 7B (1) | no data | () | () | () | no data | 1B (5) - 5A (2) |
Details of all QTL
Let’s save a huge table that will contain absolutely every QTL found as significant in this study. It will be a supplementary data for the publication.
OR=rbind(summary_QTL_DON, summary_QTL_NOT, summary_QTL_PEPI, summary_QTL_NEPIL, summary_QTL_PEPIL)All signif QTL for publication
In the publication we will talk about QTL on blups only, when available, and on raw data if not. Let’s list of these QTLs:
OR=OR[which(OR$carac%in%PUBLI) , ]
dim(OR)## [1] 72 10
OR=OR[order(OR$chromo, OR$carac), ]
dim(OR)## [1] 72 10
OR$position=round(OR$position,2)
write.table(OR, file="~/Dropbox/Publi_Fusariose/SUPPORTING_DATA/OR_detail_all_QTLs.csv", row.names=F, col.names=T, quote=F, sep=";", dec=,)Save environment
save(list=ls() , file="/Users/yan/Dropbox/Publi_Fusariose/ANALYSIS_REPRO/DATA/QTL_R_environment.R")